Coupled cluster benchmarks of water monomers and dimers extracted from DFT 
liquid water: the importance of monomer deformations 
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To understand the performance of popular density-functional theory (DFT) exchange-correlation 
(xc) functionals in simulations of liquid water, water monomers and dimers were extracted from a 
PBE simulation of liquid water and examined with coupled cluster with single and double excitations 
plus a perturbative correction for connected triples [CCSD(T)]. CCSD(T) reveals that most of 
the dimers are unbound compared to two gas phase equilibrium water monomers, largely because 
monomers within the liquid have distorted geometries. Of the three xc functionals tested, PBE and 
BLYP systematically underestimate the cost of the monomer deformations and consequently predict 
too large dissociation energies between monomers within the dimers. This is in marked contrast 
to how these functionals perform for an equilibrium water dimer and other small water clusters in 
the gas phase, which only have moderately deformed monomers. PBEO reproduces the CCSD(T) 
monomer deformation energies very well and consequently the dimer dissociation energies much more 
accurately than PBE and BLYP. Although this study is limited to water monomers and dimers, the 
results reported here may provide an explanation for the overstructured radial distribution functions 
routinely observed in BLYP and PBE simulations of liquid water and are of relevance to water in 
other phases and to other associated molecular liquids. 



I. Introduction 

Density-functional theory (DFT) has been widely used 
to study liquid water. However, how well DFT with pop- 
ular exchange-correlation (xc) functionals such as PBE 
[l| and BLYP 0, [1[ performs in describing the structural 
and dynamic properties of liquid water is a matter of 
more than a little contention. The debates, which are 
numerous, have hinged on issues such as the radial dis- 
tribution functions (RDFs) (in particular the 0-0 and 
0-H RDFs), diffusion coefficient, and average number of 
hydrogen bonds (HBs). 

It is now clear that most standard DFT molecular dy- 
namics (MD) simulations with PBE and BLYP predict an 
overstructured RDF compared to experiment. By over- 
structured, we mainly mean that the first peak in the 
0-0 RDF (referred to as ffo-o) ^ s hig ner than experi- 
ment. Consequently the computed diffusion coefficient 
is too small and the average number of HBs too large. 
Extended discussions on the magnitude and origin of the 
overstructuring can be found in e.g. Refs. 
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In brief, some of the relevant factors include: (i) The in- 
trinsic error associated with a given xc functional (in- 
cluding an improper account of van der Waals forces 



1£ i 26|); (ii) The omission of quantum nuclear effects 
H EE IS S EH; and (iii) The simulation protocol, 
with relevant factors in this regard being: (a) number of 
water molecules in the simulation cell |16j ] ; (b) the de nsity 
of the water within the cell @,H,|{J; (c) basis set [Iol.l32|: 
(d) fictitious electron mass in Car-Parrinello MD simula- 



tions [HI ; and so on. Since the first DFT MD simulation 
of liquid water in 1993 [34j . important strides have been 
made to understand how each of the above factors impact 
upon the computed properties of liquid water. However, 
simultaneously addressing all issues that could account 
for the difference between the experimental and theoret- 
ical RDFs and diffusion coefficients is not practicable, 
not to mention the uncertainties that are present in the 
experimental data itself [H, Hi| |37|. Therefore, it has 
become common to attempt to shed light on the perfor- 
mance of DFT xc functionals for treating water by inves- 
tigating well defined gas phase water clusters for which 
precise comparison can be made to high level quantum 
chemistry methods. This approach has been useful and 
allowed the intrinsic accuracy of many xc functionals to 
be precisely established M, MM\MM\MMMM, 
51, \M\ElL HI, IHJ information that may be of 
relevance to liquid water. 
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With few exceptions [47], |48|, |49|, |5fJ, |51| , previous gas 
phase benchmark studies of water clusters have focussed 
on exploring equilibrium or other stationary point config- 
urations of the gas phase intermolecular potential energy 
surfaces. However, in the liquid the structures of water 
clusters and even the water monomers themselves can 
be considerably different from those of gas phase clus- 
ters. For example, the distribution of intramolecular O- 
H bond lengths in the liquid ranges from ~0.75 to ^1.25 
A [H,[33. Yet in gas phase water clusters such as dimers 
to hexamers 0-H bond lengths deviate by <0.05 A from 
the equilibrium water monomer 0-H bond length of 0.96 
A Hi, [H| . Whether or not the performance of DFT 
xc functionals obtained from gas phase studies on wa- 
ter clusters holds for the 'deformed' structures present 
in the liquid (referred to throughout this article as 'de- 
formed') remains an important open question. Indeed 
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there is already evidence that the benchmark reference 
data obtained from gas phase clusters does not easily 
translate to the liquid. For example, BLYP predicts a 
dissociation energy for the equilibrium gas phase water 
dimcr that is 35 meV too small, yet at a water density of 1 
g/cm 3 it predicts a 3o-o that ^ s a bout 5%-15% too high. 
Similarly, PBE predicts the dimer dissociation energy to 
within 10 meV precision, yet yields an even greater ^o-o 
than BLYP. Related to this, MD simulations of liquid wa- 
ter have shown that the computed ffo-o can ^ e consider- 
ably reduced if the O-H bonds in the water monomers in 
the liquid are held rigid at some predefined bond length. 
Specifically, Allesch et al. found that the PBE gg-O de- 
creased by ~10% upon going from fully relaxed water 
monomers to a liquid with monomer O-H bonds fixed at 
~1 A [53]. Likewise, Leung et al. have shown through 
a careful and systematic series of simulations that the 
length of the O-H bonds for rigid water MD simulations 
directly correlates with [54j : as the intramolecular 

O-H bonds are allowed to lengthen, <7o-o increases. 

The studies with rigid water and the realization that 
water monomer and cluster structures in the liquid are 
likely to differ considerably from gas phase water clusters 
has prompted us to assess the performance of DFT 
xc functionals on water structures more representative 
of those present in the liquid. To this end we report 
herein on the accuracy of three DFT xc functionals 
for various deformed monomers and dimers taken from 
a PBE simulation of liquid water. Two of the most 
popular generalized gradient approximation (GGA) 
xc functionals for liquid water simulations (PBE and 
BLYP) and one of the most accurate hybrid functionals 
for small water clusters (PBEO [111) are assessed here. 
As a reference, coupled cluster with single and double 
excitations plus a perturbative correction for connected 
triples [CCSD(T)] is used with energies extrapolated 
to the complete basis set limit (CBS). The CCSD(T) 
reference calculations reveal that 75% of the dimers 
extracted from within the first coordination shell of 
the liquid are unbound relative to two equilibrium (gas 
phase) water monomers. This is mainly due to the large 
deformation of the monomers inside the liquid compared 
to the gas phase equilibrium monomer structure. PBE 
and BLYP consistently underestimate the cost of the 
monomer deformation, specifically, O-H bond stretching. 
As a consequence, both PBE and BLYP systematically 
overbind the deformed dimers extracted from the liquid, 
by as much as 80 and 43 meV, respectively. These errors 
are much larger than the usual errors associated with 
these xc functionals for the gas phase equilibrium dimer 
[3^ . In general, the performance of PBEO is superior 
to the two GGAs but noticeable errors are identified for 
all functionals including PBEO for the particularly long 
O-H bonds encountered at the shortest 0-0 separations. 
Although this study is restricted to monomers and 
dimers (and in a sense resembles a highly limited cluster 
expansion study of the liquid), the results reported here 
provide a possible explanation for the overstructured 
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FIG. 1: (Color online) From a PBE MD simulation of 
liquid water, water dimers are extracted (e.g. the high- 
lighted dimer in yellow). Single point energy calculations are 
then performed with CCSD(T), PBE, BLYP, and PBEO on 
the deformed dimers (-Edimer) and the constituent deformed 
monomers (Ei, E^). These energies are then used to evalu- 
ate the electronic dissociation energy of the dimers (Eqn. (4)) 
and the associated 1-body (Eqn. (2)) and 2-body (Eqn. (3)) 
energies. The deformation of the monomers compared to a 
gas phase equilibrium monomer is also quantified with (Eqn. 
(I))- 



RDFs routinely observed in BLYP and PBE simulations 
of liquid water. The significance of these results to water 
in other phases and to other associated molecular liquids 
is also briefly discussed. 

II. Methods, procedures, and definition of 
parameters 

Several methods have been employed to study the 
water clusters examined here. We now briefly describe 
some of the relevant computational details, how the 
water monomers and dimers are selected from the liquid, 
and then define the energetic and structural parameters 
used in the subsequent analysis. 

A. Liquid water 

To generate water monomer and dimer structures rep- 
resentative of those present in liquid water a Born- 
Oppenheimer molecular dynamics simulation of 32 D2O 
molecules in a periodic cubic box of length 9.8528 A 
was performed with the CPMD code [56]. The PBE xc 
functional was used along with hard pseudopotentials of 
Goedecker et al. [571 ] and an associated plane wave energy 
cut-off of 125 Ry. This simulation was run for 30 ps with 
an integration time step of 0.5 fs. A Nose-Hoover chain 
thermostat was used to maintain a target temperature of 
330 K. 

Water monomers and dimers were then extracted 
from the MD simulation. To get an uncorrelated 
sample of structures, 6-7 dimers were selected each 2 ps 
over the last 20 ps of the MD trajectory. In total 66 
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FIG. 2: (a) Distribution of the number of water dimers se- 
lected from the PBE M D simulation of liquid water (plotted as 
bars) as a function of the O-O separations (Ro-o) within the 
dimer. (b) Distribution of the number of O-H bonds for the 
water monomers selected from the liquid (plotted as bars). 
The dashed lines represent the corresponding O-O and O- 
H RDFs from the same PBE liquid water simulation. Only 
dimers within the first coordination shell of liquid water are 
considered here. 



that are widely employed in DFT simulations of liquid 
water (PBE and BLYP), and the hybrid PBEO func- 
tional, which is one of the most accurate functionals in 
predicting the absolute dissociation energies of small gas 
phase water clusters (dimer to hexamer) |26l . |38| . 

The CCSD(T) calculations were performed with the 
NWChem code [6(| with localized Gaussian basis sets. 
Specifically, aug-cc-pVnZ basis sets (n = T, Q, and 5) 
were used and the resultant energies extrapolated to the 
complete basis set limit (CBS) with the same standard 
heuristic schemes as employed by us before [HI, [3e| . 
CCSD(T)/CBS is the theoretical 'gold standard' for 
systems of the size considered here and, in the following, 
differences between a given xc functional and CCSD(T) 
are referred to as errors with that xc functional. In total 
>600 CCSD(T) calculations have been performed for 
the reference data presented in this paper. 

C. Definition of Parameters 

In order to quantitatively compare the structure of the 
molecules extracted from the liquid to a gas phase equi- 
librium monomer we define a quantity Sd, the deforma- 
tion, as, 



bonded dimers were selected (comprising 92 individual 
monomers). The criteria we followed for selecting 
dimers were that: (i) they were from within the first 
coordination shell of the 0-0 RDF, i.e., all chosen 0-0 
distances (i?o-o) are <3.4 A; and (ii) the distribution 
of all 66 Ro-o of the dimers resembles the 0-0 RDF 
for the first coordination shell. Fig. Ufa) illustrates 
that the distribution of dimers selected does indeed 
resemble the computed 0-0 RDF of our MD simulation 
reasonably well. As an independent check we note that 
the distribution in the values of the intra-molecular 
O-H bond lengths associated with all the selected water 
molecules is also in reasonably good agreement with the 
first peak of our computed O-H RDF of liquid water 
[Fig. Mb)}- 

B. Clusters 

The PBE water monomers and dimers extracted from 
the liquid water simulation were then examined with a 
few DFT xc functionals and CCSD(T). Throughout this 
study structures collected from the liquid water simula- 
tions are used for single point energy calculations with 
the various methods and no geometries are optimized, 
unless explicitly stated otherwise. The DFT calculations 
on the gas phase water monomers and dimers were per- 
formed with the GAUSSIAN03 code [H[, using large 
Dunning correlation consistent aug-cc-pV5Z basis sets 
[591 ]. We have shown before that for DFT xc function- 
als such as the ones considered here, this basis set is 
large enough to get dissociation energies within about 1 
meV/H 2 of the CBS limit [38]. Results from just three 
xc functionals will be reported, specifically two GGAs 



S d = /^(Rjv-rjv) 2 
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where N is the number of atoms, R and r denote the 
coordinate vectors of deformed and gas phase equilibrium 
monomer structures, respectively. 

Several energy terms will appear repeatedly and it is 
also useful to define them here. The one-body energy 
(-Eif,) of a water monomer is calculated as, 



E±b — E% -^equilibrium 



(2) 



where -Equilibrium is the energy of the gas phase water 
monomer at equilibrium and Ei is the energy of a de- 
formed monomer. The two-body energy (E 2 b) of a dimer 
is defined as: 



E 2 b — E, 



dimer 



— Ei — Eij 



(3) 



where -E^mer is the total energy of the dimer. The elec- 
tronic dissociation energy (D e ) of the dimers is given by, 



D P = E, 



dimer 



-2xE 



equilibrium 



(4) 



Fig. Q] schematically illustrates each of the above ener- 
getic quantities and the overall procedure used in this 
study. Since the structures considered here have been 
taken from a PBE liquid water simulation, -Ecquiiibrium is 
calculated with a PBE structure [6l|. The error con- 
ceded by DFT xc functionals (AE) in comparison to 
CCSD(T)/CBS is given as, 



AE = E, 
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FIG. 3: (Color online) (a) CCSD(T) one-body energy (E lb ) 
versus deformation (Sd) for water monomers taken from PBE 
liquid water. The horizontal and vertical dashed lines indi- 
cate the mean value of En, and the mean value of the defor- 
mation, respectively, (b) Comparison of En, of PBE, BLYP, 
and PBE0 with CCSD(T). 



where -Eccsd(t) an d -Edft are energies obtained from 
CCSD(T)/CBS and DFT, respectively. 

III. Results 
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FIG. 4: (Color online) (a) Variation in the one-body energy 
(Eib) with the O-H bond lengths of a water monomer cal- 
culated with CCSD(T), PBE, BLYP, and PBE0. The inset 
shows the differences in Eib of the three xc functionals com- 
pared to CCSD(T). (b) Errors in E lb (AE lb ) for the deformed 
monomers selected from liquid water as a function of the 
longest O-H bond of each monomer. The vertical dashed line 
indicates the gas phase equilibrium O-H bond length (0.97 A) 
of a monomer (optimized with PBE) and the horizontal solid, 
dashed, and dash-dotted lines represent the average errors of 
PBE, BLYP, and PBE0, respectively. Here, a positive error 
of Eib indicates that it is too easy to stretch O-H bonds of the 
monomers with a given xc functional compared to CCSD(T). 



Now we examine the water monomers extracted from 
the liquid, focusing on the cost to go from the gas 
phase equilibrium monomer structure to the deformed 
structures present in the liquid. Following this we 
consider the water dimers. In each case we compare 
the results of the various DFT xc functionals to the 
CCSD(T)/CBS references. 

A. Monomers 

To begin, CCSD(T) was used to establish the relative 
energies of the monomers taken from the liquid compared 
to the gas phase equilibrium monomer, i.e., CCSD(T) 
-E16 energies were computed for all 92 monomers. As can 
be seen from Fig. these are distributed in a very 

large range from ~0 to +900 meV with a mean value 
of +147 meV. Thus on average the monomers extracted 
from the liquid are 147 meV less stable than the gas 
phase equilibrium monomer, a surprisingly large energy. 
In quantifying the amount of deformation (Eqn. [1]) for 
each monomer we find, as expected, a general increase 
in E\h with the extent of deformation [Fig. OJa)]. The 
average deformation of the monomers is ~0.1 A, which 
gives us a measure of how deformed water monomers are 
in a PBE liquid water structure. 

Now we consider how the deformation energies 
computed with PBE, BLYP, and PBE0 compare to 
CCSD(T). This is shown in Fig. E^b), a parity plot of 
En, for the three xc functionals compared to CCSD(T). 
Immediately it can be seen that the performance of the 



GGAs is markedly different from the hybrid PBE0 func- 
tional. Specifically, for PBE and BLYP, Ei b is systemat- 
ically too small compared to CCSD(T). On average the 
PBE and BLYP deformation energies are 32 and 38 meV, 
respectively, smaller than that obtained from CCSD(T). 
The size of the error simply increases with the total 
CCSD(T) E u [Fig. [3tb)] and for the largest deforma- 
tions is on the order of 100 meV. Remembering that the 
monomer deformation is, of course, an endothermic pro- 
cess, smaller values of E\ b therefore indicate that it is too 
easy to deform monomers to their liquid structures with 
PBE or BLYP compared to CCSD(T). In contrast to the 
GGAs, PBE0 produces E\ b in excellent agreement with 
CCSD(T) with a mean error of only -3 meV. This small 
negative error reveals that it is marginally too expensive 
to deform the monomers to their liquid water structure 
with PBEO compared to CCSD(T). Since the only differ- 
ence between PBE and PBEO is the 25% Hartree-Fock 
(HF) exchange in the latter, we conclude that the inclu- 
sion of exact exchange remedies the large error in Eu, 
almost completely. Why this is so will be discussed in 
section IV. 

The monomers extracted from the liquid have both 
modified bond lengths and H-O-H internal angles. In or- 
der to understand in detail where the errors in Eu for the 
GGAs come from we carried out a simple series of tests 
where bond lengths and the internal angle were varied 
independently. The tests show that the main error in 
the GGAs comes from the bond stretching. For exam- 
ple Fig. H£a) shows that the symmetric stretching of the 
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TABLE I: Computed harmonic vibrational frequencies for a 
water monomer. V\ and vi are the asymmetric and symmetric 
O-H stretching modes and is the H-O-H bending mode. All 
values are in cm -1 and calculated with an aug-cc-pVTZ basis 
set. 









V3 


CCSD(T) 


3921 


3812 


1648 


PBE 


3800 


3696 


1593 


BLYP 


3756 


3655 


1596 


PBEO 


3962 


3856 


1633 



O-H bonds of a water monomer costs much less energy 
with PBE and BLYP compared to CCSD(T). The errors 
increase almost linearly with the stretching [inset Fig. 
HJa)] an d are as large as —200 meV when the O-H bonds 
are 0.16 A longer than the gas phase equilibrium bond 
length of 0.97 A. A bond stretch of 0.16 A may sound 
like a lot but monomers with O-H bonds even as long as 
1.18 A are present in our PBE MD simulation and in ex- 
periment and in ab initio path integral simulations even 
longer O-H bonds are observed [U H3] . As we saw for the 
structures taken from the liquid, PBEO is in very good 
agreement with CCSD(T) and even for the longest O-H 
bond of 1.18 A comes within 34 meV of CCSD(T). In ad- 
dition alterations of the H-O-H angle were considered but 
this makes much less of a contribution to the error in the 
xc functionals than what we find for bond stretching. For 
example, increasing (decreasing) the bond angle by 15° 
causes a maximum error of 15 meV (-8 meV) with BLYP 
and even smaller errors for the two other functionals. 

The above tests establish that an inaccurate descrip- 
tion of bond stretching is the main origin of the error 
in En, for the GGAs. Returning to the structures taken 
from the liquid we therefore plot in Fig. Effb) the Eu 
error against the length of the longest O-H bond for each 
monomer. As with the systematic deformations of the 
equilibrium monomer, the errors in En, increase almost 
linearly with the O-H bond length for the GGAs and 
for PBEO they remain very close to zero except at the 
longest distances. Thus it can be inferred that monomers 
inside liquid water are energetically too easy to stretch 
for both PBE and BLYP. We note that a careful series of 
tests taking water clusters from a BLYP MD simulation 
along with subsequent tests with CCSD(T) established 
that none of the conclusions arrived at here are altered 
if BLYP structures are used [62j • 

Before moving on to the dimers, we point out that 
the discrepancies established here between the three 
xc functionals and CCSD(T) correlate well with the 
errors in the computed harmonic vibrational frequencies 
of an isolated water monomer (Tabic [TJ) . Specifically 
for the two stretching frequencies PBE and BLYP are 
-115 to -160 citT 1 (-3-4%) softer than CCSD(T) 
(Table H}, whereas, PBEO is only -45 cm" 1 (-1%) 
harder. This one to one correspondence between error 
in harmonic vibrational frequencies and En, may also 
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FIG. 5: (a) CCSD(T) total dissociation energies (D e ) of the 
water dimers taken from the PBE liquid water simulation and 
the dissociation energy curve for a fully optimized gas phase 
dimer (Eqm.) as a function of the 0-0 distance (Ro-o)', (b) 
CCSD(T) two-body energies (-E2&) as a function of Ro-o for 
the same dimers. The dashed lines represent the mean values 
of D e and E-2b in panels (a) and (b), respectively. 



hold for other xc functionals and may therefore provide 
a cheap diagnostic to estimate in advance how reliably 
an xc functional will be at the determination of Eu- 

B. Dimers 

Now we move to the dimers extracted from the liquid, 
discussing first what CCSD(T) reveals about the stabil- 
ity of the dimers and then considering how well the three 
functionals perform. In Fig. 03(a) the CCSD(T) dissocia- 
tion energies are plotted as a function of the 0-0 distance 
within each dimcr. Also reported is the equilibrium (i.e., 
fully optimized) CCSD(T) dissociation energy curve for 
a gas phase water dimer. As expected the equilibrium 
dimer binding energy curve provides a lower bound for 
the dissociation energies of the deformed dimers, which 
at each particular value of Ro-O exhibit a range of val- 
ues reflecting the range of dimer structures in the liquid. 
More importantly, Fig. Et^a) provides an overview of the 
range of dissociation energies for water dimers found in- 
side the first coordination shell of PBE liquid water. The 
range is large: from -95 to +993 meV, with the mean 
value being +201 meV. Indeed 75% of the dissociation 
energies are positive, i.e., 75% of the dimers are unbound 
compared to two gas phase equilibrium water monomers. 
Upon decomposing the dissociation energies into the one- 
and two-body contributions we find that the average En, 
is 339 meV and the average E 2 b is -137 meV (Table ILl]) . 
Note that the average value of En, for the dimers is, 
of course, about twice En, for the monomers discussed 
above. Also note that this is a considerably larger En, 
than for the gas phase equilibrium water dimer, which 
is only 10 meV (Table [n}. The two-body energy gives 
the binding between the water molecules and 97% of the 
dimers have an attractive E 2 b [Fig. E£b)]. The average 
value of E-2b at -137 meV is somewhat smaller than the 
corresponding value for the gas phase equilibrium dimer 
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TABLE II: Mean values of the one-body, two-body, and dissociation energies of the deformed dimers (selected from the liquid) 
and the corresponding values for the gas phase equilibrium dimer with CCSD(T), PBE, BLYP, and PBEO. The differences 
between each DFT xc functional and CCSD(T)/CBS are given in parenthesis. Energies are in meV. 



deformed 


Gas phase equilibrium 




D e (AD e ) 


Ei b (AEi b ) 


E 26 (AE 26 ) 


D e (AD e ) 


Ei 6 (AEib) 


E 2b (AE 2 b) 


CCSD(T) 


201.9 


339.2 


-137.2 


-211.6 


9.7 


-221.4 


PBE 


121.9 (80.0) 


268.0 (71.2) 


-146.1 (8.8) 


-219.9 (8.3) 


3.7 (6.0) 


-223.6 (2.2) 


BLYP 


159.3 (42.6) 


254.4 (84.8) 


-95.1 (-42.1) 


-178.7 (-32.9) 


2.4 (7.3) 


-181.1 (-40.3) 


PBEO 


200.7 (1.2) 


346.1 (-6.9) 


-145.3 (8.1) 


-213.5 (1.9) 


9.7 (0.0) 


-223.2 (1.8) 
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FIG. 6: (Color online) (a) Errors in dimer dissociation en- 
ergy (AD e ) as a function of the sum of the deformation of 
the two monomers within each dimer for PBE, BLYP, and 
PBEO. (b) Errors in the two-body energy (AE 2b ) from PBE, 
BLYP, and PBEO as a function of the O-O distance (Ro-o) 
within the dimers. Horizontal solid, dashed, and dash-dotted 
lines represent the average errors of PBE, BLYP, and PBEO, 
respectively. 



of -221 meV. Since the total dissociation energy for the 
dimers is just the sum of E^ and E 2 b, it is quite obvious 
therefore that En, plays the major role in destabilizing 
the dimers. 

Coming back to the performance of the xc functionals, 
Fig. HJa) reports the error in the dissociation energies for 
each dimer. It can be seen that overall PBEO performs 
very well, yielding an average error of 1 meV. BLYP and 
PBE, on the other hand, yield quite large average errors 
of 43 and 80 meV, respectively. This behavior differs sig- 
nificantly from how these two functionals perform for the 
gas phase equilibrium water dimer, where errors of only 
-33 and 8 meV are obtained (Table IIT|) . Thus we arrive at 
a key result, the performance of the two GGA functionals 
for the deformed dimer structures is inferior to what it is 
for the gas phase equilibrium dimer, with both function- 
als substantially overbinding the dimers taken from the 
liquid. Table ITI1 reports the key quantities that allow us 
to understand these results and why they contrast to the 
equilibrium gas phase dimers. As one might anticipate 
from section III A, the key is the one-body deformation 



energy. In the gas phase equilibrium dimer the abso- 
lute value of Eib is small (10 meV) and the resultant 
errors even smaller (Table [TTJ) . Thus the performance of 
a functional for the gas phase equilibrium dimer is dom- 
inated by E2b, which is accurately described with PBE 
and PBEO (8 and 2 meV errors, respectively) and under- 
estimated by some 33 meV with BLYP. However, as we 
have seen in the structures taken from the liquid, En, is 
large [339 meV from CCSD(T)] and the associated errors 
from the GGA xc functionals become significant. Specif- 
ically, since both BLYP and PBE predict that the one- 
body deformation energy is much too small and predict 
either too weak or about right two-body energies then 
the total dissociation energies come out too large. BLYP 
proves to be more accurate than PBE simply because of 
more favorable cancelations of errors in E\b and E 2 b (Ta- 
ble [TTJ) . Since BLYP is generally considered to produce 
too weak HBs between water molecules [H, 0, [45[ 
the overbinding observed here is remarkable. The obvi- 
ous relevance of this finding to liquid water will be dis- 
cussed below. 

A characteristic feature of HBs between water 
molecules is that the covalent O-H bonds of the donor 
molecules (i?o-H d ) are elongated (63[. The elongation 
is a result of charge transfer from the acceptor water 
molecule to the O-H a* antibonding orbital of the donor 
molecule [64], HBj]. Since we find that it is too easy to 
stretch an O-H bond with the GGAs, one can anticipate 
that this will further influence the strength of the HBs 
formed and in particular E 2 b- To investigate this we 
return to the gas phase equilibrium water dimer as a 
test case and systematically stretch the O-H^ bond 
whilst keeping all other atoms fixed. Fig. [7{a) plots the 
change in E 2 b as a function of the O-H^ bond length 
with CCSD(T) and the three xc functionals. Clearly all 
methods predict that as the O-H^ bond increases so too 
does E 2 b- However, all three xc functionals predict too 
rapid an increase compared to CCSD(T). This is best 
seen by the inset in Fig[7[a) which displays the error in 
the change of E 2 b as a function of i?o-H d compared to 
CCSD(T). Likewise, E 2 b increases slightly too rapidly 
with the three xc functionals for the dimers extracted 
from the liquid [Fig. 0Jb)]; this is particularly apparent 
for PBE and PBEO. Thus in addition to it being too 
easy to stretch an O-H bond with BLYP and PBE, for 
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FIG. 7: (Color online) (a) Compared to the gas phase equilib- 
rium dimer, the change in the two-body energy (i?2i>) obtained 
from the systematic variation of the covalent O-H bond of the 
donor 'H' atom [5(Ro-H d )] keeping all other atoms fixed. The 
inset shows the difference between the DFT xc functionals and 
CCSD(T). (b) Error in the two-body (E 2b ) energy from DFT 
compared to CCSD(T) as a function of the 7?o-H d , obtained 
from the dimers from liquid water. Here positive values refer 
to a stronger two-body interaction. The dashed, dotted, and 
dash-dotted lines are quadratic fits to the PBEO, PBE, and 
BLYP data, respectively. 



all three xc functionals the magnitude of the change in 
E2b upon stretching is too great, further contributing to 
the overbinding of dimers with long O-H^ bonds. 

IV. Discussion 

It is clear from the last section that, despite the E^b 
errors for the longest 0-H^ bonds, the overall perfor- 
mance of PBEO is superior to that of the GGAs. To 
understand the origin of the difference between PBE and 
PBEO the variation in the exchange and correlation en- 
ergies was examined upon going from the gas phase equi- 
librium to the deformed water monomer structures (Ta- 
ble |TTl]). The variations in DFT exchange and correla- 
tion energies are then compared with the full HF exact 
exchange and CCSD(T) correlation. We note that the 
physical interpretation of exchange and correlation dif- 
fers from DFT (PBE) to CCSD(T) and so use the data 
reported in Table IIIII and Fig. [5] merely in the hope 
of obtaining some general qualitative insight. The ba- 
sic finding from CCSD(T) is that upon going from the 
gas phase equilibrium to the deformed monomers there 
is a gain in the (negative) correlation energy and a loss 
in exchange energy [Fig. Hfa),(b)]. Naturally, the abso- 
lute change in the exchange energy is far greater than 
that in the correlation energy. The two DFT xc func- 
tionals predict a loss in the exchange energy but in con- 
trast to CCSD(T) also a loss in the correlation energy, 
i.e., there is less correlation with PBE in the deformed 
monomers compared to the equilibrium monomer in the 
gas phase. Thus, in terms of the correlation energy, 
PBE/PBE0 predicts qualitatively different behavior from 



TABLE III: Average differences in the exchange and correla- 
tion contributions between the deformed monomers extracted 
from liquid water and the gas phase equilibrium monomer 
structure, obtained with CCSD(T), PBE, and PBEO. Note 
that the exchange contribution of CCSD(T) refers to HF exact 
exchange. Values in parenthesis are the differences between 
the two xc functionals and CCSD(T). Here positive and neg- 
ative values indicate energy loss and gain, respectively. All 
values are in meV. 





CCSD(T) 


PBE 


PBEO 


Correlation 


-65.5 


+35.2 (-99.7) 


+34.5 (-99.0) 


Exchange 


+825.6 


+681.9 (+143.7) 


+722.7 (+102.9) 


Total 


761.1 


717.1 (+44.0) 


757.2 (+3.9) 



CCSD(T) upon bond stretching. However, the missing 
correlation in the deformed monomers is compensated for 
by differences from CCSD(T) in the exchange energy. In 
PBEO, which predicts exchange energies in better agree- 
ment with HF (CCSD(T)) [Fig. Mb)}, the missing corre- 
lation is compensated by missing exchange so that overall 
the total energy changes are, as we have seen, very simi- 
lar for PBEO and CCSD(T). In PBE, however, the lack of 
correlation is not sufficient to compensate for the larger 
underprediction of the exchange energy [Fig. E^b)]. Thus 
we find that PBEO is superior to PBE simply because of a 
more favorable cancelation of the differences of exchange 
and correlation from CCSD(T) [Fig. [5(c)] . 

More generally the poor description of covalent O-H 
bond stretching observed here with PBE and BLYP is 
likely to apply to many other GGAs. For example, tests 
with RPBE, mPWLYP, and BP86 (on the 92 monomers 
taken from our MD simulation) all produce one-body en- 
ergies that are 30-40 meV smaller than CCSD(T). Simi- 
larly, as with PBEO, the hybrid functionals B3LYP and 
X3LYP predict rather accurate one-body energies, com- 
ing within 10 meV of CCSD(T). Of course the conclusion 
reached here that HF exact exchange is necessary for the 
proper description of covalent bond stretching is consis- 
tent with what has long been known in the context of 
covalent bond breaking (and transition state energies) in 
the gas phase (see e.g. Refs. (rMl^.lrMl^.l70l|). 

Several previous studies have examined how DFT 
xc functionals perform in treatin g H Bs between water 
molecules in water clusters [H [II Ell El El H El, El, 
E3, El, El, HI El|. Our study here is somewhat uncon- 
ventional in that we have examined structures extracted 
directly from a liquid water simulation instead of explor- 
ing equilibrium gas phase structures. This has revealed 
significant differences in how PBE and BLYP perform 
for the structures extracted from the liquid compared to 
the known performance of these functionals for equilib- 
rium gas phase water clusters. Thus we see that the 
behavior of these functionals for e.g. the gas phase equi- 
librium water dimer is not a good indicator for how these 
functionals perform for dimers extracted from the liquid. 
One must always exercise caution in making connections 
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i?2b for a gas phase equilibrium water dimer suggests that 
the likely errors in Eu for this non-hybrid GGA will be 
more than compensated for. In addition, our results are 
consistent with and help to explain the results from the 
rigid water MD simulations [53, |54|. First, by fixing O-H 
bonds at or close to the gas phase monomer equilibrium 
O-H bond length, the En, error is eliminated or greatly 
reduced. Second, the large En, error associated with the 
longest O-Hrf bonds is also obviated. 

Water molecules in other environments such as those 
in bulk ice or larger gas phase clusters will also possess 
deformed monomers with elongated bonds. These defor- 
mations are smaller than in liquid water but the effect 
is not negligible. For example, the average deformation 
of the monomers in a water hexamer is 0.05 A with 
PBE optimized geometries and in bulk ice Ih PBE pre- 
dicts an average deformation of —0.06 A. Based on the 
approximate relation between Eu> error and deformation 
established in Fig. [31 such deformations as encountered 
in small clusters and ice are likely to lead to errors in En, 
of -30-40 meV. 

Finally, we point out that the suggestion that too 
facile bond stretching may result in an overstructured 
liquid is likely to be of relevance to other associated 
liquids apart from water. The relevant experimental 
and theoretical RDFs of other associated liquids are not 
as well established as liquid water. However, there are 
indications of BLYP simulations yielding overstructured 
RDFs for e.g. liquid ammonia [7l|, [72| ancl methanol 
[73l |74| despite BLYP underestimating the strength of 
the corresponding gas phase dimers by —45% compared 
to CCSD(T) 0. 



FIG. 8: (Color online) Variation in (a) correlation energy 
(A£ corrc i a tion) and (b) exchange energy (A£ GXC han gc ) with Ct- 
rl bond length for a gas phase water molecule, (c) Variation 
in the difference of DFT correlation and exchange energies 
in comparison to the CCSD(T) correlation and HF exact ex- 
change energies, respectively, as a function of the monomer 
O-H bond length. Positive and negative values refer to energy 
loss and gain, respectively. 



between interaction energies of gas phase clusters and 
RDFs for the corresponding liquid phase, particularly in 
the present circumstances where we have only consid- 
ered one- and two-body terms. Nonetheless, it is plau- 
sible that the overbinding observed here for PBE and 
BLYP, which originates in E±b errors, is connected with 
the overstructuring of these functionals for liquid water. 
Indeed, because of the greater error cancelations between 
the one- and two-body energies calculated with BLYP, 
the overbinding of the dimers from the liquid is less for 
BLYP than for PBE. This may again provide an expla- 
nation for why <?o°o ^ s ^ ess m BLYP compared to PBE. 
This thinking may also explain the low value of #o-o 
reported in MD simulations with the RPBE functional 
[Tg. l2lj. A computed ~50 meV underestimation of the 



V. Summary 

In summary, from a PBE simulation of liquid water, 
monomers and bonded dimers (from the first coordi- 
nation shell of the 0-0 RDF) were extracted. With 
CCSD(T) 75% of the dimers were shown to be unbound 
compared to two gas phase equilibrium water monomers. 
This is mainly because the structures of the water 
monomers inside the liquid differ significantly from an 
equilibrium gas phase monomer. Indeed, with CCSD(T) 
we find that the average monomer extracted from 
the PBE liquid is about 150 meV less stable than an 
equilibrium gas phase water monomer. Among the three 
xc functionals tested, the two GGAs (BLYP and PBE) 
underestimate the energy cost for monomer deformation 
(i.e., Eif,) and as a consequence BLYP and PBE predict 
dissociation energies that are too large by 80 and 43 
meV, respectively, compared to CCSD(T). This is 
inferior to the performance of these functionals for the 
equilibrium water dimer and other water clusters in the 
gas phase. Overall PBEO yields much more accurate 
dimer dissociation energies, mainly because it is not 
susceptible to such large bond stretching errors as the 
GGAs are. However, PBEO is not free from deficiencies 
in treating the dimers examined here. Specifically, like 
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the two other functionals, it predicts an increasing error 
in E 2 b for the longest O-H^ bonds. Finally, we have 
discussed the possible relevance of these results to DFT 
simulations of liquid water, to water in other environ- 
ments, and to other associated liquids. In particular, we 
have suggested that the overbinding identified here may 
provide an explanation for the overstructured RDFs 
observed in BLYP and PBE simulations of liquid water. 
However, more work is required to further test this 
suggestion, with e.g. larger clusters that give access to 
higher order terms in the many-body decomposition 
and/or clusters embedded in an external electrostatic 
field that mimics the remaining water present in the 



liquid. 
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